library(DAseq)
library(Seurat)
Attaching SeuratObject
library(ggplot2)
library(patchwork)
library(stringr)
# set the python path
python2use <- "~/anaconda3/bin/python3"
my.samplename <- "Gg_ctrl_poly_int"
my.se <- readRDS("~/spinal_cord_paper/data/Gg_ctrl_poly_int_seurat_070223.rds")
DimPlot(my.se, reduction = "tsne", label = TRUE)
cell.labels <- unname(substr(my.se$orig.ident, 1, 4))
resp <- "poly"
non_resp <- "ctrl"
# Get the DA scores for all cells
da.cells <- getDAcells(
X = my.se@reductions$pca@cell.embeddings[, 1:30],
cell.labels = cell.labels,
labels.1 = resp, #responder labels (blue)
labels.2 = non_resp, # nonresponder labels (red)
k.vector = seq(50, 500, 50),
plot.embedding = my.se@reductions$tsne@cell.embeddings
)
Calculating DA score vector.
Running GLM.
Test on random labels.
Setting thresholds based on permutation
# same for UMAP
da.cells.u <- getDAcells(
X = my.se@reductions$pca@cell.embeddings[, 1:30],
cell.labels = cell.labels,
labels.1 = resp, #responder labels (blue)
labels.2 = non_resp, # nonresponder labels (red)
k.vector = seq(50, 500, 50),
plot.embedding = my.se@reductions$umap@cell.embeddings
)
Calculating DA score vector.
Running GLM.
Test on random labels.
Setting thresholds based on permutation
With a random permutation ofn the labels, a null distribution is generated. Min and max of the permutation DA set the default threshold.
# look at output
str(da.cells[1:2])
List of 2
$ cell.idx: int [1:15537] 1 2 3 4 5 6 7 8 9 10 ...
$ da.ratio: num [1:15537, 1:10] -0.4171 0.0939 -0.0272 0.0939 0.0939 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : NULL
.. ..$ : chr [1:10] "50" "100" "150" "200" ...
# plot the vis of the DA scores
p1 <- da.cells$pred.plot +
ggtitle(paste0("Responder \"", resp, "\" (blue) and Non-responder \"", non_resp,"\" (red)"))
p2 <- DimPlot(my.se,
ncol = 4,
reduction = "tsne",
pt.size = 0.1,
split.by = "orig.ident")
p3 <- DimPlot(my.se,
group.by = "orig.ident",
reduction = "tsne",
pt.size = 1)
p4 <- da.cells.u$pred.plot +
ggtitle(paste0("Responder \"", resp, "\" (blue) and Non-responder \"", non_resp,"\" (red)"))
p5 <- DimPlot(my.se,
ncol = 4,
reduction = "umap",
pt.size = 0.1,
split.by = "orig.ident")
p6 <- DimPlot(my.se,
group.by = "orig.ident",
reduction = "umap",
pt.size = 1)
p2
p3 + p1
p5
p6 + p4
# plot randomisation of DA scores
da.cells$rand.plot +
ggtitle(paste0("Responder \"", resp, "\" (blue/down) and Non-responder \"", non_resp,"\" (red/up)"))
With updateDAcells we can apply our own DA threshold to be more (or less) stringent.
# update DA cells with custom threshold
da.cells <- updateDAcells(
X = da.cells,
pred.thres = c(-0.7,0.7),
plot.embedding = my.se@reductions$tsne@cell.embeddings
)
# check output
da.cells$da.cells.plot +
ggtitle(paste0("Responder \"", resp, "\" (blue) and Non-responder \"", non_resp,"\" (red)"))
da.cells.u <- updateDAcells(
X = da.cells,
pred.thres = c(-0.7,0.7),
plot.embedding = my.se@reductions$umap@cell.embeddings
)
# check output
da.cells.u$da.cells.plot +
ggtitle(paste0("Responder \"", resp, "\" (blue) and Non-responder \"", non_resp,"\" (red)"))
# get the contiguous regions of cells
da.regions <- getDAregion(
X = my.se@reductions$pca@cell.embeddings[, 1:20],
da.cells = da.cells,
cell.labels = cell.labels,
labels.1 = resp,
labels.2 = non_resp,
resolution = 0.01,
plot.embedding = my.se@reductions$tsne@cell.embeddings,
)
Using min.cell = 50
Warning: The following arguments are not used: row.names
Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')
Warning: The following arguments are not used: row.names
Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')
Warning: k.param set larger than number of cells. Setting k.param to number of
cells - 1.
Warning: The following arguments are not used: row.names
Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')
Removing 4 DA regions with cells < 50.
str(da.regions[1:2])
List of 2
$ cell.idx : int [1:15537] 1 2 3 4 5 6 7 8 9 10 ...
$ da.region.label: num [1:15537] 0 0 0 0 0 0 0 0 0 0 ...
da.regions$da.region.plot
# get the contiguous regions of cells (UMAP)
da.regions.u <- getDAregion(
X = my.se@reductions$pca@cell.embeddings[, 1:20],
da.cells = da.cells,
cell.labels = cell.labels,
labels.1 = resp,
labels.2 = non_resp,
resolution = 0.01,
plot.embedding = my.se@reductions$umap@cell.embeddings,
)
Using min.cell = 50
Warning: The following arguments are not used: row.names
Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')
Warning: The following arguments are not used: row.names
Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')
Warning: k.param set larger than number of cells. Setting k.param to number of
cells - 1.
Warning: The following arguments are not used: row.names
Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')
Removing 4 DA regions with cells < 50.
da.regions.u$da.region.plot
# add da regions as meta data to plot them
my.se$da.regions <- da.regions$da.region.label
# cell x gene matrix
expr <- as.matrix(my.se[["RNA"]]@data)
STG.markers <- STGmarkerFinder(
X = expr,
da.regions = da.regions,
lambda = 1.5, n.runs = 5, return.model = T,
python.use = python2use
)
# wilcoxon rank sum (MAST does not work)
Seurat.markers <-SeuratMarkerFinder(
object = my.se,
features = my.se[["integrated"]]@var.features,
da.slot = "da.regions",
da.regions.to.run = 1:(length(table(da.regions$da.region.label))-1),
min.pct = 0.25,
logfc.threshold = 0.5,
return.thresh = 0.05,
assay = "RNA"
)
Save the output:
saveRDS(da.cells, paste0("~/spinal_cord_paper/output/",my.samplename,"_DA_cells.rds"))
saveRDS(da.regions, paste0("~/spinal_cord_paper/output/",my.samplename,"_DA_regions.rds"))
saveRDS(da.cells.u, paste0("~/spinal_cord_paper/output/",my.samplename,"_DA_cells_umap.rds"))
saveRDS(da.regions.u, paste0("~/spinal_cord_paper/output/",my.samplename,"_DA_regions_umap.rds"))
saveRDS(STG.markers, paste0("~/spinal_cord_paper/output/",my.samplename,"_STG_markers.rds"))
saveRDS(Seurat.markers, paste0("~/spinal_cord_paper/output/",my.samplename,"_Seurat_markers.rds"))
pdf(paste0("~/spinal_cord_paper/figures/",my.samplename,"_DAseq.pdf"), width = 10, height = 10)
DimPlot(
my.se,
group.by = "da.regions",
reduction = "tsne",
cols = c(
"grey",
rainbow(
length(table(da.regions$da.region.label))-1
)
),
label = TRUE,
repel = TRUE,
label.size = 6
)+
ggtitle(paste0(my.samplename," DA regions"))+
theme(plot.title = element_text(hjust = 0.5))
DimPlot(
my.se,
group.by = "da.regions",
reduction = "umap",
cols = c(
"grey",
rainbow(
length(table(da.regions$da.region.label))-1
)
),
label = TRUE,
repel = TRUE,
label.size = 6
)+
ggtitle(paste0(my.samplename," DA regions"))+
theme(plot.title = element_text(hjust = 0.5))
p2
p1
p5
p4
da.cells$rand.plot
da.cells$da.cells.plot +
ggtitle(paste0("Responder \"", resp, "\" (blue) and Non-responder \"", non_resp,"\" (red)"))
da.cells.u$da.cells.plot +
ggtitle(paste0("Responder \"", resp, "\" (blue) and Non-responder \"", non_resp,"\" (red)"))
dev.off()
png
2
# Date and time of Rendering
Sys.time()
[1] "2023-02-27 12:12:45 CET"
sessionInfo()
R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)
Matrix products: default
BLAS/LAPACK: /scicore/soft/apps/OpenBLAS/0.3.1-GCC-7.3.0-2.30/lib/libopenblas_sandybridgep-r0.3.1.so
locale:
[1] en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] stringr_1.4.0 patchwork_1.1.1 ggplot2_3.3.3 SeuratObject_4.0.2
[5] Seurat_4.0.5 DAseq_1.0.0
loaded via a namespace (and not attached):
[1] plyr_1.8.6 igraph_1.2.6 lazyeval_0.2.2
[4] sp_1.4-5 splines_4.1.0 listenv_0.8.0
[7] scattermore_0.7 digest_0.6.27 foreach_1.5.1
[10] htmltools_0.5.1.1 fansi_0.5.0 magrittr_2.0.1
[13] tensor_1.5 cluster_2.1.2 ROCR_1.0-11
[16] limma_3.48.0 recipes_0.1.16 globals_0.16.2
[19] gower_0.2.2 matrixStats_0.58.0 spatstat.sparse_3.0-0
[22] colorspace_2.0-1 ggrepel_0.9.1 xfun_0.34
[25] dplyr_1.0.10 jsonlite_1.7.2 spatstat.data_3.0-0
[28] survival_3.2-11 zoo_1.8-9 iterators_1.0.13
[31] glue_1.6.2 polyclip_1.10-0 gtable_0.3.0
[34] ipred_0.9-11 leiden_0.3.9 future.apply_1.7.0
[37] shape_1.4.6 abind_1.4-5 scales_1.1.1
[40] DBI_1.1.1 miniUI_0.1.1.1 Rcpp_1.0.7
[43] viridisLite_0.4.0 xtable_1.8-4 reticulate_1.22
[46] spatstat.core_2.1-2 proxy_0.4-25 stats4_4.1.0
[49] lava_1.6.9 prodlim_2019.11.13 glmnet_4.1-1
[52] htmlwidgets_1.5.3 httr_1.4.2 RColorBrewer_1.1-2
[55] ellipsis_0.3.2 ica_1.0-2 farver_2.1.0
[58] pkgconfig_2.0.3 nnet_7.3-16 sass_0.4.0
[61] uwot_0.1.10 deldir_1.0-6 here_1.0.1
[64] utf8_1.2.1 caret_6.0-88 tidyselect_1.2.0
[67] labeling_0.4.2 rlang_1.0.6 reshape2_1.4.4
[70] later_1.2.0 munsell_0.5.0 tools_4.1.0
[73] cli_3.4.1 generics_0.1.3 ggridges_0.5.3
[76] evaluate_0.20 fastmap_1.1.0 yaml_2.2.1
[79] goftest_1.2-2 ModelMetrics_1.2.2.2 knitr_1.41
[82] fitdistrplus_1.1-6 purrr_0.3.4 RANN_2.6.1
[85] pbapply_1.4-3 future_1.30.0 nlme_3.1-152
[88] mime_0.10 compiler_4.1.0 plotly_4.10.0
[91] png_0.1-7 e1071_1.7-7 spatstat.utils_3.0-1
[94] tibble_3.1.8 bslib_0.2.5.1 stringi_1.6.2
[97] highr_0.9 lattice_0.20-44 Matrix_1.3-3
[100] vctrs_0.5.1 pillar_1.8.1 lifecycle_1.0.3
[103] spatstat.geom_3.0-3 lmtest_0.9-38 jquerylib_0.1.4
[106] RcppAnnoy_0.0.19 data.table_1.14.0 cowplot_1.1.1
[109] irlba_2.3.3 httpuv_1.6.1 R6_2.5.0
[112] promises_1.2.0.1 KernSmooth_2.23-20 gridExtra_2.3
[115] parallelly_1.33.0 codetools_0.2-18 MASS_7.3-54
[118] assertthat_0.2.1 rprojroot_2.0.2 withr_2.4.2
[121] sctransform_0.3.3 mgcv_1.8-35 parallel_4.1.0
[124] grid_4.1.0 rpart_4.1-15 timeDate_3043.102
[127] tidyr_1.1.3 class_7.3-19 rmarkdown_2.17
[130] Rtsne_0.15 pROC_1.17.0.1 shiny_1.6.0
[133] lubridate_1.7.10